
# find the index of the first 1 in the row vector v
indexOfFirstOne := proc(V)
local i, v;
    i := 1;
    for v in V do
        if v = 1 then
            return i;
        end if;
        i := i + 1;
    od;
    i
end proc;

# perform column reordering to achieve systematic form
adjustColumn := proc(H)
local rownum, colnum, ind, A, i;
local pms;
    pms := [];
    (rownum, colnum) := LinearAlgebra:-Dimension(H);
    A := copy(H);
    for i from 1 to rownum do
        ind := indexOfFirstOne(A[i]);
        if (ind <> i) then
            print("switching columns of H:");
            print(i);
            print(ind);
            pms := [[i, ind], op(pms)];
            LinearAlgebra:-ColumnOperation(A, [i, ind], inplace=true);
        end if;
    od;
    [pms, A]
end proc;


# This function puts parity check matrix H into the form [ I | A ]
convertH := proc(H)
    adjustColumn(Gaussjord(H) mod 2);
end proc;

# This function computes the generator matrix of a parity check
# matrix.
computeG := proc(H)
local H2, P, PT, rownum, colnum, A, i, j, k, pm, pms, pair, R;
# n - k = rownum, n = colnum
    (rownum, colnum) := LinearAlgebra:-Dimension(H);
    k := colnum - rownum;
    pair := convertH(H);
    pms := pair[1];
    H2 := pair[2];
    PT := LinearAlgebra:-SubMatrix(H2, [1..rownum], [colnum-k+1..colnum]);
    P := LinearAlgebra:-Transpose(PT);
    R := ArrayTools:-Concatenate(2, P, Matrix(k, shape=identity));
    # perform inverse column permutation to switch the columns back
    for pm in pms do
        LinearAlgebra:-ColumnOperation(R, pm, inplace=true);
    od;
    R
end proc;






# a testing case from Wikipedia:
# A := Matrix([[1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 0, 1], [1, 0, 0, 1, 1, 0]]);
# H := convertH(A);
# G := computeG(H);

# creates the first submatrix H_1. The matrix has n/k rows. Each row
# has k 1's, descending.
getDescendingMatrix := proc(n, k)
local M, rownum, colnum, i, j;
    colnum := n;
    rownum := iquo(n, k);
    M := Matrix(rownum, colnum, 0);
    for i from 1 to rownum do
        for j from (i-1)*k + 1 to i*k do
            M(i, j) := 1;
        od;
    od;
    M
end proc;

with(combinat, randperm):

# computes the matrix whose columns are generated from the column
# permutations of the input matrix M.
permuteColumns := proc(M, n)
local cols, R, i;
    cols := randperm([LinearAlgebra:-Column(M, [1..n])]);
    # R := cols[1];
    # for i from 2 to n do
    #     R := ArrayTools:-Concatenate(2, R, cols[i]);
    # od;
    Matrix(cols)
end proc;



# generate a random number r s.t. 1 <= r <= n, and r < low and r >
# high
randomButNotInRange := proc(n, S)
local r;
    r := (rand() mod n) + 1;
    while member(r, S) do
        #print(S);
        r := (rand() mod n) + 1;
    od;
    #print(r);
    r
end proc;


getLow := proc(r, k)
local a;
    a := iquo(r, k, 'x');
    if x = 0 then return r - k + 1 end if;
    a * k + 1
end proc;

# computes the matrix whose columns are generated from the column
# permutations of the input matrix M.
permuteColumnsRst := proc(M, n, k)
local cols, R, i, r, low, high, S, S2, j;
    cols := [LinearAlgebra:-Column(M, [1..n])];
    #for the first column, randomly choose one.
    r := rand(1..n)();
    R := cols[r];
    low := getLow(r, k);
    high := low + k - 1;
    S := {seq(low..high)};
#    print("init set S");
#    print(S);
    S2 := {r};
    for i from 2 to n do
        r := randomButNotInRange(n, S);
        R := ArrayTools:-Concatenate(2, R, cols[r]);
        low := getLow(r, k);
        high := low + k - 1;
#        print("previous S");
#        print(S);
        S := `union`(S, {seq(low..high)});
#        print("new S");
#        print(S);
        S2 := `union`(S2, {r});
        if (i mod k) = 0  then
            S := copy(S2);
        end if;
    od;
    Matrix(R)
end proc;

# find the smallest cycle with length 4 in the parity check matrix.
containCycle4 := proc(U, V, n)
local i;
    for i from 1 to (n-1) do
        if U[i]*V[i]*U[i+1]*V[i+1] = 1 then
            return true;
        end if;
    od;
    false
end proc;

overlapMoreThanOne := proc(H, P, r)
local rownum, colnum, i, v, j;
    (rownum, colnum) := LinearAlgebra:-Dimension(H);
    for i from 1 to r do
        v := P[i];
        for j from 1 to rownum do
            if H[j] . v > 1 then
                return true;
            end if;
        od;
    od;
    false
end proc;

# Gallager's (n, j, k)-code construction where
# n: code length
# j: column weight
# k: row weight
# the function return a pair
# H: the original Gallager's parity check matrix.
# G: the generator matrix generated from the transformed H.
# FIXME:  currently buggy becuase length-4 cycle can not be avoided
# fully.
GallagerCodes := proc(n, j, k)
    local H, P, topMat, i, a, b;
    # compute the top descending matrix.
    topMat := getDescendingMatrix(n, k);
    H := Matrix(topMat);
    for i from 1 to (j-1) do
        P := permuteColumns(topMat, n);
        # if a length 4 cycle is found, redo the column permutation
        # matrix generation.
        #while containCycle4(LinearAlgebra:-Row(H, -1), LinearAlgebra:-Row(P, 1), n) do #or
        while overlapMoreThanOne(H, P, n/k) do
            P := permuteColumns(topMat, n);
        od;
        # concat to the matrix on the top
        H := ArrayTools:-Concatenate(1, H, P);
    od;
    [H, computeG(H)]
end proc;

#HG := GallagerCodes(20, 3, 4);


#=========================== McKay's code construction ===================
# generate a set containing t random numbers.
# the numbers have to be >= 1 and <= k
genRandomNumberSet := proc(k, t)
    local S, r, i;
    # seed the random number generator with system clock.
    #Seed := randomize();
    S := {};
    for i from 1 to t do
        r := (rand() mod k) + 1;
        while member(r, S) do
            r := (rand() mod k) + 1;
        od;
        S := {op(S), r};
    od;
    S
end proc;

# generate a random column vector of weight t
genOneColumn := proc(k, t)
    local V, S, i;
    V := Vector([seq(0, i=1..k)]);
    S := genRandomNumberSet(k, t);
    for i in S do
        V[i] := 1;
    od;
    V
end proc;

# check overlap constraints to make sure there's no two columns have
# overlap more than 1. This is to avoid length-4 cycle.
has4Cycle := proc(H, V)
local rownum, colnum, i, Hi;
    (rownum, colnum) := LinearAlgebra:-Dimension(H);
    for i from 1 to colnum do
        Hi := LinearAlgebra:-Column(H, i);
        if Hi . V > 1 then
            return true;
        end if;
    od;
    false
end proc;

# compute the row weights of a parity-check matrix
# for examining the property of the parity-check matrices.
rowWeights := proc(H)
    local L, i, n, k;
    (k, n) := LinearAlgebra:-Dimension(H);
    L := [];
    for i from 1 to k do
        L := [op(L), add(H[i][j], j=1..n)];
    od;
    L
end proc;

# compute the column weight of a parity-check matrix via MacKay's
# construction.
columnWeight := proc(H)
    local n, k, w, i;
    (k, n) := LinearAlgebra:-Dimension(H);
    w := 0;
    for i from 1 to k do
        w := w + H[i][1];
    od;
    w
end proc;

# MacKay's code construction 1A.
# avoid 2-column overlapping
# the parity check matrix has dimension k * n and column weight t.
# FIXME:
# 1. need to keep on generating until a linearly independent matrix is
# generated.
# 2. need to make row weight as uniform as possible.
# 3. need to avoid row weight 1.
MacKay1A := proc(n, k, t)
local R, V, i, Seed;
    # seed the random number generator with system clock.
    Seed := randomize();
    R := Matrix(genOneColumn(k, t));
    for i from 2 to n do
        V := genOneColumn(k, t);
        while has4Cycle(R, V) do
            V := genOneColumn(k, t);
        od;
        #print("generate column:");
        #print(i);
        R := ArrayTools:-Concatenate(2, R, V);
    od;
    R  #[R, computeG(R)]
end proc;

#HG := MacKay1A(500, 300, 3);

# encode a message msg with generator matrix G
encodeLDPC := proc(msg, G)
    convert(LinearAlgebra[Modular]:-Multiply(2, Vector[row](msg), G), list)
end proc;
